function fjac = VecFun_Jacobian(f,param,x)

tol     =   eps.^(1/3);
h       =   tol.*max(abs(x),1);
xh1     =   x+h; 
xh0     =   x-h;
h       =   xh1-xh0;
f1      =   f(xh1,param);
f0      =   f(xh0,param);
fjac    =   (f1-f0)./h;